# title: wang-ho2010 estimation algorithm
# author: StableGenius
# date: 2020-01-23
# description: 复现wang-ho2010提出的估计算法
# references:
## - Wang HJ, Ho CW. Estimating fixed-effect panel stochastic frontier models by model transformation[J]. Journal of Econometrics, 2010, 157(2):286-296.

# import packages
library(Formula)
library(plm)
library(tidyverse)
library(frontier)

# drop intercept in x and z
drop_intercept <- function(formula){
  if (is_formula(formula)){
    formula <- Formula(formula)
  }
  update(formula, .~.-1|.-1)
}

# within transformation
within_trans <- function(data, index_panel, panel_index){
  tmp_panel <- pdata.frame(cbind(index_panel, data), 
                           index = panel_index)
  
  result <- matrix(NA, nrow = nrow(index_panel), ncol = ncol(tmp_panel)-2)
  for (index in 1:ncol(result)) {
    result[, index] <- Within(tmp_panel[[index+2]])
  }
  
  result
}

# Wang-ho2010 estimation
sfa_wh <- function(formula, data){
  formula <- drop_intercept(formula)

  # check panel data
  if (!"pdata.frame" %in% class(data)){
    stop("data set is not panel data!")
  }
  
  if (!pdim(data)$balanced){
    stop("unbalanced panel data!")
  }
  
  # dimension of panel data
  dimN <- pdim(data)$nT$n
  dimT <- pdim(data)$nT$T
  
  panel_index <- names(attr(data, "index"))
  index_panel <- data[, panel_index]
  
  # get original variables and do within transformation
  y <- model.response(model.frame(formula, data)) %>%
    within_trans(index_panel, panel_index)
  X <- model.matrix(formula, data, rhs=1) %>%
    within_trans(index_panel, panel_index)
  K <- ncol(X)
  Z <- model.matrix(formula, data, rhs=2)%>%
    within_trans(index_panel, panel_index)
  Kz <- ncol(Z)
  
  # parameters initialization
  fe_model <- plm(formula = formula(terms(formula, lhs=1, rhs=1)), 
                  data = data, 
                  effect = "individual", 
                  model = "within")
  beta_0 <- coef(fe_model)
  
  sfa_sigmaSq <- var(residuals(fe_model))
  sigma_u2 <- 0.1 * sfa_sigmaSq
  sigma_v2 <- sfa_sigmaSq - sigma_u2
  delta_0 <- rep(0, Kz)
  
  parameters <- c(beta_0, delta_0, sigma_u2, sigma_v2)
  names(parameters) <- c(names(beta_0), 
                         attr(terms(formula, rhs=2), "term.labels"),
                         "sigmauSq", "sigmavSq")
  # log likelihood
  logL <- function(param){
    p_beta <- param[1:K]
    p_delta <- param[(K+1) : (K+1+Kz)]
    p_sigmau2 <- param[length(param)-1]
    p_sigmav2 <- param[length(param)]
    
    if (p_sigmav2<0 | p_sigmau2<0)
      return(NA)
    
    ep <- (y - X %*% p_beta) %>%
      as.vector()
    h <- Z %*% p_delta %>%
      as.vector()
    
    lT <- rep(1, dimT)
    Pi <- (diag(dimT) - lT %*% t(lT) /dimT) * p_sigmav2
    Pi_inv <- MASS::ginv(Pi)
    
    mu <- 0
    result <- 0
    
    for (index in 1:dimN) {
      ep_i <- ep[(dimT * (index-1) + 1) : (index * dimT)]
      h_i <- h[(dimT * (index-1) + 1) : (index * dimT)]

      mu_star <- (mu/p_sigmau2-t(ep_i) %*% Pi_inv %*% h_i) /
        (t(h_i) %*% Pi_inv %*% h_i + 1/p_sigmau2)
      sigma2_star <- 1/(t(h_i) %*% Pi_inv %*% h_i + 1/p_sigmau2)
      sigma_star <- sqrt(sigma2_star)
      
      
      result <- result + 
        -(dimT-1)/2*log(2*pi) +
        -(dimT-1)/2 * log(p_sigmav2) +
        -0.5*t(ep_i) %*% Pi_inv %*% ep_i +
        0.5*(mu_star^2/sigma2_star - mu^2/p_sigmau2) + 
        log(sigma_star * pnorm(mu_star/sigma_star)) +
        -log(sqrt(p_sigmau2)*pnorm(mu/sqrt(p_sigmau2)))
    }
    
    result
  }
  
  maxLik(logLik = logL, 
         start = parameters,
         method = "NM")
}